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The collective excitations in ensembles of dissipative, laser driven ultracold atoms exhibit crystal¬ 
like patterns, a many-body effect of the Rydberg blockade mechanism. These crystalline structure 
are revealed in experiment from a post-selection of configurations with fixed numbers of excitations. 
Here, we show that these sub-ensemble can be well represented by ensembles of effective particles 
that interact via logarithmic pair potentials. This allows one to study the emergent patterns with a 
small number of effective particles to determine the phases of Rydberg crystals and to systematically 
study contributions from iV-body terms. 
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Current experiments with laser driven, trapped ultra¬ 
cold Rydberg atoms m open a viable route to study 
driven dissipative many-body systems in the strongly cor¬ 
related regime [zhu]. In these experiments atoms are 
prepared in the ground state and are collectively driven 
by a laser to a Rydberg state. Due to dipole-dipole or 
van-der Waals interactions, excitations of adjacent parti¬ 
cles are suppressed, the socalled Rydberg blockade mech¬ 
anism [Fig. [5], In ongoing experiments, the position 
of individual Rydberg excitations can be determined in 
single-shot single-particle resolved measurements [Tj [a 
typical configuration is illustrated in Fig. [5ja)] . These 
excitation configurations are then post-selected into sub¬ 
ensembles according to the number of excitations N, re¬ 
sulting in a conditional density p(x\N), where x repre¬ 
sents the positions of all N excitations [see Fig[5|[b)]. It is 
these sub-ensembles that exhibit the remarkable crystal¬ 
like density patterns known as Rydberg crystals. While 
the short time coherent regime has been studied experi¬ 
mentally, we will focus here on the steady state including 
dissipation mm- As the blockade radius may be large 
compared to the the lattice spacing, the number of atoms 
N a can be several orders of magnitude larger than the 
number of excitations N. 

In this Letter we present an analytic effective parti¬ 
cle model for excitations in dissipative Rydberg crystals. 
This model is derived from the conditional density, which 
can be written as a sum of ./V-body terms 

p(x|IV) = exp[d>(x)] = (1) 

= exp [$ 0 + ^$ 1 (x i ) + ^$ 2 (a;j,^) + 

i ij 

+ X ^3{Xi,Xj,X k ) + ...]• 

ijk 

Here, the first term $0 is a global shift, $ 1 ( 2 ^) repre¬ 
sents a term that only depends on single particle posi¬ 
tions, $ 2 (xi,Xj) is a pair correlation term, etc. Below 
we will present an analytic expression for <k 2 with a log¬ 
arithmic distance dependence. We refer to this term as 



FIG. 1: (a) The system of interest consists of ground-state 
atoms (red dots) which are laser excited to Rydberg states 
(green dots) in an optical lattice (see e.g. Ref. [I]). Exci¬ 
tations within a blockade radius rb are suppressed, (b) Con¬ 
ditional density p(x\N = 7) from the logarithmic potential 
model, (c) Excitation probability (a rr ) (see main text) as a 
function of the distance r in units of the excitation maximum 
r m = (A/C 6 ) 1/6 for A = 40. For resonant (A = 0, black) 
and red-detuned driving A < 0 (red) excitations are sup¬ 
pressed at short distances. When driving with a blue-detuned 
laser A > 0 (blue) the excitation probability is peaked at a 
resonant distance r m . (d) Logarithmic potential Eq. § as 
a function of the distance r between two effective particles. 
In the case A = 0,—40,—30,—20,—10 (black and red), the 
effective particles are purely repulsive. For blue detunings 
A = 10, 20, 30, 40 (blue) the excitation resonance translates 
into a potential energy minimum. (Other parameters in units 
of O are: C 6 = l,r K = 0.1, T z = 1). 


logarithmic pair interaction in analogy to particle inter¬ 
actions in statistical mechanics [see Fig. [fjjd)]. We show 
that the steady state of Rydberg crystals p(x|iV) is well 
described by the interaction <h 2 for all N. In particular, 
for van-der Waals interactions it is the dominant term. 
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We find that higher order contributions show a spatial 
pattern which indicate that these terms soften the crys¬ 
tals (see Fig. [4]). Results are compared with the steady 
state solution of the quantum master equation [Eq. [3] in¬ 
cluding quantum correlations for small systems and with 
a Monte Carlo scheme from Ref. [T^ for large system 
sizes, respectively. 

Model - We consider a laser driven system of ultracold 
trapped atoms described by the Hamiltonian 


N a 


h = J2 


N u 

2 K 9r 


n 


( a gr + a gr ) + Vi. 


o-Wo-0') ( 2 ) 

ijU rr o rr . 


Here, N a is the number of atoms, A is the detuning, 
(Txy = \x) (y\ are the matrix elements coupling states 
x,y £ {r,g}, where g is the ground state and r is the 
Rydberg-state, evaluated for particle i and the bare Rabi 
frequency is H. The interaction is Vij(r ) = Ce/r ®- (van- 
der Waals) or Vij(r) = C^/r^j (dipole-dipole) with the 
distance between atom i and j. In the following energies 
are given in units of f 2 and distances in units of a o, the 
lattice spacing between atoms. We are interested in the 
steady state in the open dissipative regime described by 
the master equation 

^ = -i[H,p]+f^{£ip + Cip) ( 3 ) 

i=o 

with = ^T r [2a J rg p<T J rg -{a 3 rr , p}\ and C{p = T z [(a J rr - 
a ig)p{ a ir ~ aJ gg) ~ p\ (details see [H]). In a typical ex¬ 
periment, the decay rate is small T r -C Q. Dephasing 
results from various mechanisms such as finite linewidth 
of the excitation laser and particle motion which can be 
of the order of F, ss fi. The steady state of the quan¬ 
tum master equation can be calculated numerically exact 
for small numbers of atoms (typically N a < 20). An ex¬ 
ample of p(x) is shown in Fig. [2] (a) for N a = 5. The 
conditional density is then deduced from tracing out all 
configurations with excitations p(x.\N) = Tr[p<5(A r — N)], 
where N = ’Yf"' (orr) and S is the Dirac delta function. 

Logarithmic Pair Interaction - In the following we sys¬ 
tematically study the expansion Eq. (|T]) . As the simplest 
approximation, the constant term <f>o and all higher or¬ 
der terms $ 3 , $ 4 ,... are neglected. In this case, the only 
remaining term is the pair potential <I> 2 . Note, that an ex¬ 
ternal trapping potential can be added using the $1 term. 
Further we assume that the term <f>2 (a:*, Xj ) = i>(|:Ej — Xj\) 
is a function of the distance only. In this approximation 
the steady state probability of finding a configuration x 
in the ensemble with N excitations is 


-FV(x) 
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Z N 


N 

exp[J^ $ 2 (a.’i - 
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*j)]‘ 


( 4 ) 




FIG. 2: (a) Steady-sate density p(x) of excitations from the 
quantum master equation [Eq. [ 3 ] (squares) compared to the 
results of the Monte Carlo scheme from Ref. 116] (circles) and 
the logarithmic potentials (solid line). Here, N a = 5, 12 = 1, 
A = 0, ryrv = 50 and rb = 1 with a van-der Waals direct in¬ 
teraction V(r) oc 1/r 6 . (b) Pair correlation function of excita¬ 
tions from logarithmic potentials (solid line) and Monte Carlo 
(circles). For both, van-der Waals interaction V(r) oc 1/r 6 
(black) and dipole-dipole interaction V(r) oc 1/r 3 (red), the 
results are in excellent agreement. For comparison, a crystal 
that results from the direct interaction 1/r 6 or 1/r 3 , respec¬ 
tively, has completely different features (dotted). The param¬ 
eters are chosen as in Ref. [16]: f2 = 1, A = 0, r r = 0.112, 
r z = 12. The blockade radius was set to rj = 20 oq. 


Here, the probability is normalized with the the partition 
sum Zjg. In this approximation, the probability of states 
Eq. Q has the form of a canonical distribution. In the 
following we derive a analytic expression for the potential 
$ 2 - 

For two particles, it has been shown m that 
in the large time limit t r r ,r 2 the steady 

state probability of two atoms at distance r is 
(,Prr) = D 2 /{2H 2 + ^( 7 l g + [A + V(r)} 2 )}. Here, 

7 rg = \Tr + 2 r z with F r the decay rate, F z the 
dephasing and f2 the Rabi frequency, 

where N r is the number of atoms within a blockade 
radius mm- In the case of van der Waals interactions, 
V(r) = Ce/r^ if i and j are excited and V{r) = 0 
otherwise illustrated in Fig. [5](a) together with the 
excitation probability as a function of the distance 
( p rr ). The length-scale of the a superatom is given 
by the blockade distance which is defined for van-der 
Waals interactions as rb = {Cq/{ 2VL)^/T~YJlrg) 1 ^ T61 . 
Applying the reversible work theorem |19| to connect 
the pair correlation with a pair interaction we find 
$ 2 (rij) = — ln[t 7 rr (r)] + hi[oy r (oo)] which can be written 
as 


® 2 {rij) = In 


(^) 2 + 2 A^ 
_ 12 _ 

+ 7r 2 9 + a 2 


( 5 ) 


Here, r^ = 27 — Xj is the distance between the effec¬ 
tive particle i and j. Examples of this potential with 
detunings in the range A = [—40,40] are depicted in Fig. 
§d). This particle model allows one to interpret the pat¬ 
terns of post-selected Rydberg excitations as the result of 
A-particles with a logarithmic pair interaction. In par- 
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FIG. 3: (a) Distribution of excitations in a ID chain of N a = 
200 atoms, A = 0, C 6 = lO.Ofl, R. = O.Olfi, R = 0.R1 The 
average number of excitations is N = 4.4 and configurations 
with N = 3, 4, 5, 6 are observed, (b-e) The density-density 
pair correlation function <72 (r, N) (dots) is almost identical to 
the result of the logarithmic potential model (solid) for each 
sub-set with N particles. 


ticular by interpreting $2 = V e g/(ksT), where T is the 
temperature, ks is the Boltzmann constant, and V e g the 
effective potential energy of the particles, the system fol¬ 
lows a canonical ensemble. For a comparison of thermal 
fluctuations and fluctuations in the microscopic picture 
see m ■ Note, that in this interpretation the tempera¬ 
ture is not a tunable parameter. We note, that for blue 
detuning, the first two orders of a series expansion at the 
minimum are similar to a Lennard-Jones potential. 

Numerical Results - For small system sizes Eq. |3| 
can be solved numerically and compared to the logarith¬ 
mic potential model. For large system sizes this is not 
feasible and we compare our results with results from 
a Monte Carlo scheme introduced in Ref. m which is 
shortly summarize here for completeness. In this scheme 
each configuration is represented N a atoms where N are 
excited. In each Monte Carlo step, transitions between 
configurations from N excitations to either N— 1 or IV +1 
excitations are possible with the following protocol: a 
random atom i is chosen and, given that it is in the 
ground state, it is excited with an acceptance probability 
Pace = ff 2 /( 2 S 1 2 + rg + (Ei <3 V( r ij) + A)“]). Ac¬ 
cordingly, if the chosen atom is excited it is de-excited 


with probability 1 — P aC c- Note, that the logarithmic 
potential model acts on sub-sets of constant number of 
excitations, and the sequence of de-excitation and exci¬ 
tation is translated to a dynamics of an excitation due to 
pair interactions [20]. A useful measure to compare con¬ 
figurations is the pair correlation function of excitations 

9 2 {r) = [ ri “ r j])>’ ( 6 ) 

where angular brackets (•) denote the ensemble average 
in the the canonical ensemble and the average over all 
configurations. 

The comparison for small system sizes as depicted in 
Fig. [2ja) shows a good agreement between the results 
of the quantum master equation, the Monte Carlo proce¬ 
dure from Ref. [16j and the logarithmic potential model. 
The density p(x) in the logarithmic potential model is 
calculated from the p(x) = c n p{x\n), where the sum 
goes from n = 1 to n = 5 and the coefficients c n is are 
taken from an independent Monte Carlo simulation. 

Fig. Kb) depicts the results of the pair correlation 
functions from Monte Carlo calculations in comparison 
to the logarithmic potential model for a large ID sys¬ 
tem. For the chosen parameters, the average number of 
excitations is N = 20. We find that the spatial excitation 
statistics from the dynamics of particles with logarithmic 
pair interaction is almost identical to the results from 
the Monte Carlo simulation. For comparison, the crystal 
structure from particles with a direct van-der-Waals or 
dipole-dipole interaction exhibits a completely different 
pattern. 

The logarithmic potential model does not just repro¬ 
duce the pattern from an average over all possible excita¬ 
tion numbers p{x) but also for each sub-set of ensembles 
p(x\N). Figs. [3]depict the results for all numbers of exci¬ 
tations that are possible for a given blockade radius. The 
distribution of excitations from the Monte Carlo simula¬ 
tion [Fig. Ka)] is of Gaussian shape and peaked around 
IV « 4. To compare the pair correlation function for each 
sub-set of constant number of excitations we consider 

9 2 {r,N) = -(X^(r- [ri ~ rj]))jv, (7) 

i¥=j 

where the (A)n = {A6(N — N)) denotes averages with a 
constraint on the number of excitations N. The results 
are shown in panels Figs. §b) -(e) ranging from N = 3 to 
N = 6. For small numbers of excitations the system is in 
a liquid-like state while for larger densities the structures 
are crystal like. These results from the Monte Carlo sim¬ 
ulation are in good agreement with the results from the 
canonical ensemble of logarithmic pair interactions with 
the corresponding numbers of particles. 

As a two-dimensional example we study a disc-like 
setup as used in Ref. [Tj. In this setup, atoms are 
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FIG. 4: Comparison of the density of states after the rotation 
protocol from the full master equation with post selection of 
configurations with N = 6 excitations (a) and from logarith¬ 
mic potentials with N = 6 particles (b). Note, that the shape 
is the result of a mixture of a hexagon and a pentagon from 
configurations with all 6 excitations at the border and con¬ 
figurations with 5 excitations on the border an an excitation 
in the center, (c) Comparison of the steady state probability 
of excited atoms in the N = 6 manifold as a function of x 
for fixed y. The results from the logarithmic potential (solid) 
and master equation (dashed) are in good agreement, (d) 
The difference in density between the two models allows one 
to estimate the effects of three body (as well as higher order) 
terms. Here, the iV-body contributions are of the the order 
of 10~ 4 (approximately 1%) and the pattern indicates that 
these contributions destabilize crystallinity. 


five good agreement as shown in Fig. , 0 b). Comparing 
the densities quantitatively [see Fig. [4j)c)] shows that for 
van-der Waals interactions the deviation is of the order 
of 1%. These deviations are, however, not uniform in the 
system. The spatial distribution [Fig. §d)] shows that 
three-body terms destabilize the crystalline pattern. 

In conclusion, we have presented a framework that al¬ 
lows one to treat a system described by excitations and 
de-excitations as a series of dynamical many-body sys¬ 
tems with logarithmic pair interactions. This allows one 
to predict Rydberg crystal patterns with well-established 
tools from statistical mechanics.A practical advantage is, 
that the number of variables is reduced from N a to N 
variables which is considerable smaller for a large block¬ 
ade radius. The expansion of the density into potentials 
with TV-body terms also opens a way to study three body 
contributions. As an outlook, the mapping and the sta¬ 
tistical mechanics description of Rydberg crystals may be 
interesting as an extension to previous hard rod models of 
Rydberg excitations ED da. Further, the mapping from 
the excitations model to a interaction model may also 
serve as the starting point for the reverse process, where 
pair interactions are designed with a particular purpose 
(e.g. formation of a particular structure, anisotropic in¬ 
teraction) which can be mapped back on a excitation pro¬ 
tocol. Finally, the multi-particle terms which are small in 
the case of van-der Waals interactions may become dom¬ 
inant in systems with long range interactions and may 
open a new way to design multi-particle interactions. 
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loaded into a deep optical lattice. Then, atoms with a 
distance from the origin larger than ro are removed (see 
Fig. 0 a)]. The parameters are chosen such that there 
are approximately N a = 20 x 20 sites and the blockade 
distance rb « ro- This allows for Rydberg crystals with 
N < 8 excitations. The system is rotationally invariant 
and therefore a rotational invariant averaging protocol 
was introduced in Ref. [Tj. In this protocol, described in 
more detail in |20| , each configuration is reorientated to 
eliminate the rotational invariance. The density of exci¬ 
tations after this rotation is shown for N = 6 in Fig. |4|a) 
from Monte Carlo simulations and Fig. 0 b) from the 
logarithmic potential model. For N = 6 two symmetries 
are possible: (i) either all 6 excitations are positioned at 
the surface and the crystal is of hexagonal shape, or (ii) 
there is 1 excitation in the center is surrounded by 5 exci¬ 
tations. In the latter case the crystal has a 5-fold rotation 
symmetry. This mixture of two crystalline structures is 
clearly seen in Figs. [4|a) and also reproduced with the 
logarithmic potentials. The results are also in quantita- 
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Supplementary Materials 
Fluctuations 

The fluctuations around the equilibrium position of 
particles in the logarithmic potential model correspond 
to de-excitations of an atom succeeded by an immediate 
excitation of an other atom in the close vicinity (illus¬ 
trated in Fig. [5]). In the following we compare the dis¬ 
tribution of the location of these re-excitations with the 
thermodynamic fluctuations in the particle model. As an 
illustration we introduce the following simple model and 
analytically calculate the expected fluctuations in both 
models, the excitation model and the particle model. In 
the steady state the average distance between excitations 
is the blockade radius 



FIG. 6: The distribution of the position of a re-excited atoms 
with two neighbors at the steady state distance compared to 
the fluctuations of a particle in the logarithmic potential of 
two neighboring particles. 



FIG. 5: Illustration of the mapping: In the microscopic pic¬ 
ture, the steady state is the result of atoms that are excited 
and de-excited. Consider the following sequence: In a typical 
steady state configuration 1 (top left) with N = 7 excitations 
one Rydberg atom is de-excited (process Wt m ). In a second 
step, another atom is excited at a different position (process 
Wnm)- Assuming 1 was a steady state configuration, the prob¬ 
ability is large that the atom is excited close to the position 
of the previously de-excited atom. In the logarithmic poten¬ 
tial model, we consider only sub-ensembles of constant N, in 
this example N = 7. In this interpretation, the sequence of 
events is a movement of an effective particle following a pair 
interaction. 



We consider an excited atom at position x = 0 with two 
neighbors at x = tb and x = — rs, assuming only nearest 
neighbor interaction in a ID system. If the atom at the 
origin gets de-excited it can only be re-excited in the close 
vicinity of the original position. The distribution of the 
exact position of re-excitation is then 


Pre{ x ) — 


1 


n 2 


Qre 2fl 2 + —[7r S + i( x -r B ) 6 + (x+r B )e) 2 } 

(9) 
is the 


where Q re = V 


-+ 


^6 

0 : + T Ti) 


) 2 ] 


normalization factor. The Lindemann width of this dis¬ 
tribution is then 


Ax = 


\J■f-r B Pr e( X ) 2 
Tb 


( 10 ) 


In the particle picture, this corresponds to an atom at 
position x = 0 with thermal fluctuations in the logarith¬ 
mic potential from two particles fixed at x = —rs and 
x = tb which is 


V(x)=V e f f (x-r B ) + V e ff(x + r B ), (11) 

where V e ff = from Eq.(5) from the main 

text. The equilibrium probability of the particle is 


Ppot — n e 

Wpot 


1 

k R T 


( 12 ) 


V(a) 


Here, the normalization factor is Q po t = e 

The distributions Eq. (£ 121 and Eq. (SjsiJ for typical 
parameters (A = 0,7 r <? = 0.1D, T = 0.1H, Cq = 100H) 
are in almost perfect agreement. The Lindemann pa¬ 
rameters of these distributions are Aa: = 0.568 and 
Ai = 0.57 for the Monte Carlo simulation and particle 
picture, respectively. 


Rotation Protocol for the Spherical Geometry 

The crystalline pattern in the spherical geometry in 
Figs 4 (a) and (b) is the result of a rotation protocol to 
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allow comparison with Ref. [1]. We use the following 
method: For each excitation we calculate the distance to 
the center s = \/x 1 + y 2 and distinguish between par¬ 
ticles at the border s >= 5 and particles in the center 
s < 5. The number of excitations on the border then 
serves as the symmetry the crystal is compared to (e.g if 
6 particles are on the border Nb = 6, the crystal is com¬ 
pared to a hexagon). The parameter for each excitation 
i we calculate 

*n b =Y , eNB6i > ( 13 ) 

i 

where 6i is the angle of the vector from the origin to ex¬ 
citation i with respect to the y axis. The form of if) is a 
generalization of the hexatic order parameter x/jq which 
is used to classify the orientation of hexagonal plaquettes 
in a triangular lattice. The parameter ipN B is the relative 
orientation of a polygon with Nb side with respect to the 
polygon spanned by the Rydberg excitations. Then the 
system is rotated to find the minimum of tpN B ■ This ro¬ 
tated configuration is then added to the histogram shown 
in Figs. 4 (a) and (b) in the main text. Note, that for a 


given number of excitations several symmetries can oc¬ 
cur. For example, for N = 6 the crystal can either consist 
of 6 excitations on the border (hexagonal) or 5 excitations 
on the border (pentagonal) and 1 excitation in the cen¬ 
ter. This is clearly seen in Fig. (4) as the pattern is a 
mixture of a hexagonal and a pentagonal crystal. 


Patterns of various A'-cnsernblcs 

The excitation patterns in Rydberg ensembles are re¬ 
vealed from a post-selection of configurations with N ex¬ 
citations. This can be achieved by solving the master 
equation [Eq. (3) from the main text] and store configu¬ 
rations according to the number of excitations. The re¬ 
sulting patterns for TV = 5,6, and 7 excitations is shown 
in Fig. [7] (left column). The logarithmic potential model 
requires an individual simulation for each number of ex¬ 
citation. The potential is the same in all 3 cases. The 
results, shown in Fig. [ 7 ] (right column) are in good agree¬ 
ment with the original model for all densities. 
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FIG. 7: Series of density plots for various numbers of ex¬ 
citations. In the master equation approach (left), the full 
ensemble is split up in a post-selection step into the various 
IV-manifolds. In the logarithmic potential approach (right) 
each N corresponds to a system with N particles and loga¬ 
rithmic pair interactions. 



































